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(57)Ab5tract 

The disdosure involves the use of a library of modulated wavelet- 
packets which are effective in providing both precise frequency localiza- 
tion and space localization. An aspect of the disclosure involves feature 
extraction by determination of the correlations of a library of waveforms 
with the signal being processed, while maintaining orthogonality of the 
set of waveforms selected (i.e. a selected advantageous basis). In a dis- 
closed embodiment, a method is provided for encoding and decoding an 
input signal, comprising the following steps (80): applying combinations 
of dilations and translations of a wavelet to the input signal to obtain 
processed values (1610); computing the information costs of the pro- 
cessed values (1620); selecting, as encoded signals, an orthogonal group 
of processed values, the selection being dependent on the computed in- 
formation costs (1630); and decoding the encoded signals to obtain an 
output signal (1640). The wavelet preferably has a plurality of vanishing 
monfents. 


.1610 


GENERATE PROCESSED VALUES 
BY CORRELATING THE UBRARY OF 
WAVELET PACKETS WITH 
INPUT SIGNAL SAMPLES. 
( RGURE 17 ) 


COMPUTE INFORMATION COSTS 
OF PROCESSED VALUES. 
( RGURE 18 ) 


SELECT ADVANTAGEOUS 
ORTHOGONAL BASIS. 
( RGURE 19 ) 


1620 


1630 


1640 


COMPILE ENCODED FRAMES 
FOR STORAGE AND/OR 
TRANSMISSION. 
( RGURE 20 ) 
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Description 

METHOD AND APPARATUS FOR ENCODING AND 
DECODING USING WAVELET -PACKETS 

FIELD OF THE INVENTION 
This invention relates to a method and apparatus for 
encoding and decoding signals which may represent any continuous 
or discrete values . 

BACKGROUND OF THE INVENTION 

It is well established that various types of signals can be 
efficiently encoded and subsequently decoded in a manner which 
substantially reduces the size of the information required (e.g. 
number of bits, bandwidth, or memory) without undue or noticeable 
degradation of the decoded signal. Examples are the tj^pes of 
audio and video bandwidth compression schemes that are currently 
in widespread use. 

In signal analysis, it is often useful to recognize the 
appearance of characteristic frequencies, but this knowledge 
generally has to be coupled with the location of the time (or 
space) interval giving rise to the frequency. Such questions can 
usually be tackled by use of the windowed Fourier transform, with 
different size windows corresponding to the scale of the 
transient feature. This analysis can be achieved by correlating 
the signal to all windowed exponentials and checking for large 
correlations. Unfortunately, due to lack of independence, 
information obtained can be redundant and inefficient for feature 
reconstruction purposes . 

is among the objects of the present invention to provide 
an im* *<.ved encoding and decoding method and apparatus which 
overcf- ;53 limitations of prior art techniques and provides 
improved and more efficient operation. 

SUM^aARY OF THE INVENTION 
The present invention involves, in part, the use of a 
library of modulated wavelet-packets which are ef f ctive in 
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providing both precise frequency localization and spac 
localization. An aspect of the invention involves feature 
extraction by determination of the correlations of a library of 
waveforms with the signal being processed, while maintaining 
orthogonality of the set of waveforms selected (i.e. a selected 
advantageous basis). 

In accordance with an embodiment of the invention, a method 
is provided for encoding and decoding an input signal, 
comprising the following steps: applying combinations of 
dilations and translations of a wavelet to the input signal to 
obtain processed values; computing the information costs of the 
processed values; selecting, as encoded signals, an orthogonal 
group of processed values, the selection being dependent on the 
computed infoxnaation costs; and decoding the encoded signals to 
obtain an output signal. As used herein, wavelets are zero mean 
value orthogonal basis functions which are non-zero over a 
limited extent and are used to transform an operator by their 
application to the operator in a finite number of scales 
(dilations) and positions (translations) to obtain transform 
coefficients. [In the computational context, very small 
non-zero values may be treated as zero if they are known not to 
affect the desired accuracy of the solution to a problem. ] A 
single averaging wavelet of unity mean is permitted. Reference 
can be made, for example, to: A. Haar, Zur Theorie der 
Orthogonalen Functionsysteme, Math Annal. 69 (1910); K.G. 
Beauchamp, Walsh Functions And Their Applications, Academic 
Press (1975); I. Daubechies, Orthonormal Bases of Compactly 
Supported Wavelets, Comm. Pure Appl. Math XLl (1988). 

In a preferred embodiment of the invention, the wavelet has 
a plurality of vanishing moments. In this embodiment, the step 
'of applying combinations of dilations and translations of the 
wavelet to the input signal to obtain processed values coir^rises 
correlating said combinations of dilations and translations of 
the wavelet with the input signal . The combinations of 
dilations and translations of the wavelet are designated as 
wavelet-packets, and in a disclosed embodiment the step of 
applying wavelet-packets to the input signal to obtain process d 


wo 91/18361 


PCr/US9I/03504 


3 

values includes: generating a tr e of proc ssed values, th 
tree having successive levels obtained by applying to the input 
signal, for a given level, wavelet-packets which are 
combinations of the wavelet-packets applied at a previous level. 
Also in a disclosed embodiment, the steps of computing 
information costs and selecting an orthogonal group of processed 
values includes performing said computing at a number of 
different levels of said tree, and performing said selecting 
from among the different levels of the tree to obtain an 
orthogonal group having a minimal information cost (the "best 
basis"). Also in this embodiment, the step of selecting an 
orthogonal group of processed values includes generating encoded 
signals which represent said processed values in conjunction 
with their respective locations in said tree. 

Further features and advantages of the invention will 
become more readily apparent from the following detailed 
description when taken in conjunction with the accompanying 
drawings * 

BRIEF DESCRIPTION OF THE DRAWINGS 

Fig. 1 is a block diagram of an apparatus in accordance with 
an embodiment of the invention, and which can be used to practice 
the method of the invention. 

Fig. 2 is a diagram illustrating Haar functions and 
combinations of dilations and translations of such functions. 

Fig. 3 illustrates a tree of nodes. 

Fig.s 4A, 4B and 4C illustrate examples of possible 
orthogonal bases. 

Fig. 5 illustrates a wavelet having two vanishing moments. 

Fig.s 6-13 illustrates examples of wavelet-packets. 

Fig. 14 illustrates an example of how information cost can 
be computed. 

Fig.s 15A and 153 illustrate a procedure for reconstruction 
in accordance with an embodiment of the invention. 

Fig. 16 shows a flow diagram which is suitable for 
controlling the encoder processor to implement an embodiment of 
the encoding apparatus and m thod in accordance with the 
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invention . 

Fig. 17 is a flow diagram of a routine for generating 
processed values from the sampled signal using a wavelet-packet 
basis. 

Fig. 18 is a flow diagram of a routine for computing 
information cost. 

Fig. 19 is a routine for selecting an advantageous or best 
orthogonal basis. 

Fig. 20 is a flow diagram of a routine for generating output 
encoded words. 

Fig. 21 is a flow diagram of the decoder routine for 
processing frames of words and reconstructing the orthogonal 
basis indicated by the words of a frame. 

Fig. 22 is a further portion of the routine for 
reconstruction in the decoder. 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 
Referring to Fig. 1, there is shown a block diagram of an 
apparatus in accordance with an embodiment of the invention for 
encoding and decoding an input signal which can be any continuous 
or discrete signal or sequence of numbers representing values in 
one or more dimensions (e.g. audio, still or moving pictures, 
atmospheric measurement data, etc.) ^nd which, for purposes of 
illustration, can be considered as an audio signal x(t) . At the 
encoder 100 the signal is coupled to an analog-to-digital 
converter 102, which produces signal samples x^, x^, x^..., a 
sequence of which can be envisioned as a vector x. The digital 
samples are coupled to an encoder processor 105 which, when 
programmed in the manner to be described, can be used to 
implement an embodiment of the invention and to practice an 
embodiment of the method of the invention. The processor 105 may 
be any suitable processor, for example an electronic digital or 
analog processor or microprocessor. It will be understood that 
any general purpose or special purpose processor, or other 
machine or circuitry that can perform the computations described 
herein, electronically, optically, or by other means, can be 
utilized. The processor 105, which for purposes of the 
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particular d scribed embodiments hereof can be c nsidered as th 

processor or CPU of a general purpose electronic digital 

computer, such as a SUN-3/50 Computer sold by Sun Microsystems, 

will typically include memories 125, clock and timing circuitry 

130, input/output functions 135 and display functions 140, which 

may all be of conventional types. 

With the processor appropriately programmed, as described 

hereinbelow, a compressed output signal x*" is produced which is 

an encoded version of the input signal, but which requires less 

bandwidth. In the illustration of Fig. 1, the encoded signals x*" 

are shown as being coupled to a transmitter 140 for transmission 

over a communications medium (air, cable, fiber optical link, 

microwave link, etc.) 150 to a receiver 160. The encoded signals 

are also illustrated as being coupled to a storage medium 142, 

which may alternatively be part of t ie processor subsystem 105, 

and are also illustrated as being manipulated such as by 

multiplication by a sparse matrix M*^, as described in the U»S. 

Patent Application Serial No. 525,974. [See also Appendix III.] 

The matrix can obtain be obtained using the wavelet -packet 

best basis selection hereof (see also Appendix V). The signal 

itself may, of course, also be in the form of a matrix (i.e., a 

collection of vectors) • The stored and/or manipulated signals 

can be decoded by the same processor subsystem 105. (suitably 

programmed, as will be described) or other decoding means. 

In the illustrated embodiment, another processor 175, which 

is shown as being similar to the processor 105, also includes 

memories 225, clock and timing circuitry 230, input/output 

functions 235, and display functions 240, which may again be of 

conventional types. Processor 175 is employed, when suitably 

programmed as described, to decode the received encoded signal 

t • • 

x*", and to produce an output digital signal x^^ , x^ / X3 . . . , (or 
vector x') which is a representation of the input digital signal, 
and which can be converted, such as by digital-to-analog 
converter 195, to obtain an analog representation x'(t) of the 
original input signal x(t) . As will become understood, the 
accuracy of the representation will depend upon the encoding 
process and the degree of bandwidth compression. 
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Bef re describing the underlying theory of the invention, it 
is noted that reference can be made to Appendices I-V, appended 


foundations and further approaches. 

A well known wavelet basis, which has a single vanishing 
moment, as defined hereinbelow, is the Haar basis [see, for 
example, A. Haar, Zur Theorie Der Orthogonalen 

Funktionensysteme, Math Annal. 69, 1910, and Appendix I of the 
aboveref erenced U.S. Patent Application Serial No. 525,974]. 
Consider the Haar basis as applied to a simplified case of eight 
samples x^, x^. • • x^. For uniform amplitudes, and ignoring 
normalizing coefficients (which are multiples of 1//2 for Haar 
wavelets), a set of waveforms can be developed from combinations 
of Haar wavelets, as illustrated in Fig. 2 and in accordance with 
the following relationships: 


hereto, for supplemental description of the theoretical 



(1) 



(2) 


sds^ = ds^ 
dds^ = ds^ 


SSSj = ss^ 
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ssd^ = sd^ + sdj (3) 
dsd^ = sd^ - sdj 
sdd^ = dd^ + ddj 
ddd^ = dd^ - dd^ 

The first group of relationships, (1) [the top eight waveforms in 
Fig. 2], are Haar functions, and are orthogonal. The last group 
of relationships, (3) [the bottom eight waveforms in Fig. 2], are 
the first eight of the well known Walsh functions. As is known 
in the art, the Walsh functions are orthogonal and complete, and 
can be advantageously used to transfonn and subsequently 
back- trans form certain types of signals to achieve, inter alia, 
signal compression. It can be observed, however, that the set of 
the entire twenty- four functions of Fig. 2 is not orthogonal, 
which follows from the fact that some of the functions are 
derived from combinations of other functions. For example, sd^ 
is the sum of d^ and d^, and is not orthogonal to d^ or to d^. 

The sums and differences of relationships (l)-(3) are 
arranged in a tree of "nodes" in Fig. 3. Four "levels" are 
shown, level 0 being the sample data, and levels 1, 2 and 3 
respectively corresponding to the groups of relationships (1), 
(2) and (3) above. The boxes ("nodes") at each level contain 
respective sum and difference terms, and are connected by 
branches to the nodes from which they are derived. It is seen 
that level 1 has two nodes (labelled, from left to right, 1 and 
2), level 2 has four nodes (labelled, from left to right, 1-4), 
and level 3 has eight nodes (labelled, from left to right, 1-8). 
It follows that a level k would have 2^ nodes. The "positions" 
of the functions (or members) within a node are also numbered, 
from left to right, as illustrated in node 1 of level 1 (only). 
The nodes of level 1 each have four positions, the nodes of the 
' level 2 each have two positions, and the nodes of level 3 each 
have one position. It is seen that each "parent" node has two 
"child" nodes, the members of the children being derived from 
those of the parent. Thus, for example, node 1 of level 2 and 
node 2 of level 2 are both children of node 1 of level 1. This 
follows from the relationships (2) set forth above. In 
particular, the members of node 1, level 2 (ss^ and ss^) are 
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derived from siims of members of node 1, level 1, and the members 
of node 2, level 2 (ds^ and ds^) are derived from differences of 
members, of node 1, level 1. 

In accordance with an aspect of the present invention, a 
complete set of functions (that is, a set or "basis" which 
permits reconstruction of the original signal samples) is 
obtained from the tree, permitting selection of nodes from any 
level. The selection is made in a manner which minimizes the 
information cost of the basis; i.e., the selected basis can be 
represented with a minimum bandwidth requirement or minimum 
number of bits for a given quality of information conveyed. 
Orthogonality of the selected basis is maintained by following 
the irule that no ancestor (parent, grandparent, etc.) of a 
selected node is used. [Conversely, no descendant (child, 
grandchild, etc.) of a selected node is used.] For the 
simplified tree of Fig. 3, there are twenty five possible 
orthogonal basis selections. Using shaded boxes to indicate the 
nodes selected for a given basis, four examples of possible 
orthogonal bases are shown in Figs. 4A-4D. If desired a basis 
which is the best level basis could alternatively be determined 
and used. 

The Haar wavelet system, and wavelet-packets derived 
therefrom, has been used so far for ease of explanation. In 
accordance with an aspect of the present invention, advantageous 
wavelet-packets are generated using wavelets having a plurality 
of vanishing moments and, preferably, several vanishing moments. 
For description of the types of wavelets from which these 
wavelet -packets can be derived, reference can be made to I. 
Daubechies, Orthonormal Bases of Compactly Supported Wavelets, 
Comm. Pure, Applied Math, XLl, 1988; Y. Meyer Principe 
d' Incertitude, Bases Hilbertiennes et Alg&bres d ' Op&rateurs , 
S^minaire Bourbaki, 1985-86, 662, Ast&risque (Socidte' 
Math^matique de France); S. Mallat, Review of Multifrequency 
Channel Decomposition of Images and Wavelet Models, Technical 
Report 412, Robotics Report 178, NYU (1988), and the 
abovereferenced U.S. Patent Application Serial No. 525,974. 

The wavelet illustrated in Fig. 5 has two vanishing 
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moments. As used herein, the number of vanishing moments, for a 
wavelet 7(x) is determined by the highest integer n for which 

/^(x) x^'dx = 0 

where 0<k5n-l, and this is known as the vanishing moments 
condition. Using this convention, the Haar wavelet has 1 
vanishing moment, and the wavelet of Fig. 5 has 2 vanishing 
moments. [It is generally advantageous to utilize wavelets 
having as many vanishing moments as is practical, it being 
understood that the computational burden increases as the number 
of vanishing moments (and coefficients) increases. Accordingly, 
a trade-off exists which will generally lead to use of a moderate 
number of vanishing moments.] The wavelet of Fig. 5 has 
coefficients as follows (see e.g. Daubechies, supra): 




(1 + 

/3)/4/2 



(3 + 

/3)/4/2 



(3 - 

/3)/4i/'2 

\ 


(1 - 

/3)/4/2 

Si 








33 




94 





The vanishing moments condition, written in terms of 
defining coefficients, would be the following: 

L 

E g^i*" = 0 
i=l 

where 

05k5n-l and 
and L is the number of coefficients. 

The procedure for applying this wavelet is similar to that for 
the Haar wavelet, but groups of four elements are utilized and 
are multiplied by the h coefficients and the g coefficients to 
obtain the respective terms of opposing polarity. To obtain 
wavelet-packets from the wavelets, the previously illustrated 
levels of sum and diff rence terms are obtained using the h and g 


wo 91/18361 


PCr/US91/03S04 


10 

coefficients, respectively. The "sum" correlation terms for the 
first level are computed as follows: 

s^ = hjX^ + h^Xj + hjXj + h^x^ 


s, = h^Xj + h,x^ + hjX, + h^Xg 

(4) 


S3 = h^Xj + h,Xg + h3X^ + h^Xg 


The "difference" correlation terms for the first level are 
computed as follows: 

(5) 


The four sets of second level sum and difference correlation 
terms can then be computed from the first level values as 
follows: 

ss = h s + h s + h s + h s 
00^ "i**! "2 2 3 3 * 4 4 


ss, = h^s, + h,S, + hjS, + h^Sg 

(6) 


SS3 = h^Sg + h,Sg + hjS, + h^s^ 


= Vlc/2-1 + Vlc/2 + Vk/2*1 + Vk/2*2 


ds, = g^s^ + g,s^ + 9383 + g^s^ 
ds, = g^s, + g,s^ + g^3^ + g^s^ 
'is, = g^s^ + g,s^ + gjS, + g,s^ 


= giSk/2-1 + 92^^/2 + 93^^/2*1 + 94Sk/2« 


(7) 
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Sd, = h,d, + h,d, + h3d3 + h,d, 
Sd, = h,d3 + h,d, + h3d, + hA 
Sd, = + h,d, + h3d, + h,d3 

(8) 

dd, = g^d, + g^d, + g3d3 + g,d, 
= gA + g,d, + gjd, + g^ 
dd, = g^dj + g^ + gjd, + g^d3 

(9) 

• 

, and so on. Extra values at end positions can be handled by 
"wrap-around", truncation, or other known means of handling end 
conditions. It will be understood that the procedure described 
in conjunction with relationships (4) -(9) is operative to 
successively correlate the signal samples with the 
wavelet-packets for each successive level. If desired, the 
correlations can be implemented by generating wavelet-packets a 
priori (using the indicated coefficients, for this example), and 
then individually correlating wavelet-packets with the signal 
using analog (e.g. electronic or optical means), digital, or any 
suitable technique. If desired, a special purpose network could 
be used to perform the correlations. 

In terms of the diagram of Fig. 3, the sums (4) and the 
differences (5) would occupy nodes 1 and 2, respectively, of 
level 1, and the sums (6), differences (7), sums (8) and 
differences (9) would occupy the nodes 1, 2, 3 and 4, 
respectively, of level 2. 

It will be understood that in many practical situations the 
number of samples considered in a frame or window (level 0) will 
be larger than 8, and the tree will also be larger than those 
shown here for ease of illustration. Fig.s 6-10 show the first 
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five wavelet-packets synth sized for sampl length 1024, using a 
wavelet with six coefficients (six h*s and six g's), and Fig.s 
11-13 illustrate three of the higher frec[uency wavelet-packets. 

In accordance with an aspect of the invention, the basis is 
selected in a manner which mininiizes information cost. There are 
various suitable ways in which information cost can be computed. 
Fig. 14 shows a parent node and its two children nodes. As a 
measure of infozrmation, one can compute the number of 
correlation values in the parent node that exceed the particular 
threshold and the total number N of correlation values in the 

c 

child nodes that exceed the particular threshold. As represented 
in Fig. 14, if is less than or equal to N^, the parent node 
will be preferred, whereas if is greater than N^, the children 
nodes will be preferred. As higher level comparisons are made 
(with the ancestors of the parent) the selection may be 
supplanted by an ancestor. 

Another measure of information cost that can be used is the 
entropy cost function that is well known in information theory, 
and is threshold independent (see Appendix IV} . Suitable 
weighting of coefficients can also be used. For example, if it 
is known or computed that certain values should be filtered, 
emphasized, or ignored, weighting can be appropriately applied 
for this purpose. 

Fig. ISA illustrates a procedure for reconstruction which 
can be utilized at the decoder processor. The shaded boxes 
indicate the nodes defining the orthogonal basis that was 
selected at the encoder and is to be decoded. The arrows 
illustrate the reconstruction paths, and the cross-hatched boxes 
indicate reconstructed nodes, the last (level 0) reconstruction 
providing the reconstructed decoder output information. In 
particular, node 1, level 2 is reconstructed from node 1, level 
3 and node 2, level 3. Node 1, level 1 is then reconstructed 
from node 1, level 2 and node 2, level 2, and so on. 

Fig. 15B shows children nodes containing sy^, sy2,...sy^ and 
dy^, dy^, . . .dy^ being mapped into their parent node to obtain the 

reconstructed y^, y^, y^^. If the four coefficients h^, h^, 

and h^ (with the corresponding g^, g^, g^ and g^) were used for 


wo 91/18361 


PCT/US91/03504 


13 

encoding, the decoding relationships will be as follows: 
For V odd 

(10) 


For V even 

= hjsy^ + h.sy,, + g^dy^ + g^dy^ 
y^ = h^sy^ + h^sy^ + g^dy^ + g^dy^ 

(11) 


y2K = ^^Vn + ^Syn-l + 92<iyn + 94^yn-l 

The values mapped into the parent mode are accumulated and, as in 
the encoder, extra values at the end positions (e.g. at y^ above) 
c .r» be handled by "wrap-around", truncation, or other known means 
of handling end conditions. 

Referring to Fig. 16, there is shown a flow diagram which, 
when taken together with the further flow diagrams referred to 
therein, is suitable for controlling the processor to implement 
an embodiment of the encoding apparatus and method in accordance 
with the invention. The block 1610 represents the generating of 
processed values by correlating wavelet-packets with the input 
signal samples. There are various ways in which this can be 
achieved, the routine of Fig. 17 describing an implementation of 
the present embodiment. The block 1620 represents the routine, 
described further in conjunction with Fig. 18, for computing the 
* information costs of the processed values . As described further 
hereinbelow, there are various ways of computing the measure or 
cost of information contained in the processed values. In an 
illustrated embodiment, a thresholding procedure is utilized, and 
information cost is determined by the number of values which 
exceed a particular threshold. The block 1630 represents the 
routine of Fig. 19 for s lection of an advantageous orthogonal 
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basis from among the processed values, the s lection of the basis 
being dependent on the computed information costs . The block 
1640 represents compiling encoded frames from the selected 
processed values which constitute the basis, such as for 
subsequent recovery after processing, storage, and/or 
transmission. This routine is described in conjunction with 
Fig* 20. 

Referring to Fig. 17, there is shown a flow diagram of a 
routine for generating the processed values from the sampled 
signal, or signal portion, using a wavelet-packet basis. The 
block 1710 represents the reading in of the selected coefficients 
h^/g^. The block 1720 is then entered, this block representing 
the initializing of a level index at 0, the initializing of a 
node index at 1, and the initializing of a position index at 1. 
The scunple data, considered as level 0 is then read in, as 
represented by the block 1725. The sample data may consist, for 
example, of 256 sequential samples of an acoustical signal to be 
compressed and transmitted. The level index is then 
incremented, as represented by block 1730, and the first level 
processed values are computed and stored in accordance with the 
relationships (4) and (5) set forth above (block 1735 and loops 
1748 and 1760). For example, for the first position of the first 
node of level 1, s^ will be computed*. If the wavelet employed is 
representable by a filter having four coefficients, as in the 
example above, s^^ will be computed as the sum of h^x^^, h^x^, hjXj, 
h^x^. If a wavelet of more vanishing moments is used, more 
coefficients will be employed. In general, it will be preferable 
to utilize a wavelet having several coefficients, greater than 
four, the above exaunples being set forth for ease of 
illustration . 

In loop 1748, inquiry is made (diamond 1740) as to whether 
the last position of the current node has been reached* If not, 
the position index is incremented (block 1745), and the block 
1735 is re-entered for computation of the next processed value 
of the current node and level. The loop 1748 is then continued 
until all processed values have been computed for the current 
node, whereupon inquiry is made (diamond 1750) as to whether the 
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last nod of the current 1 vel has be n reached. If not, the 
node index is incremented (block 1755) and the loop 1760 is 
continued until the processed values have been computed for all 
nodes of the current level. For the first level, there will be 
only two nodes, with the values thereof being computed in 
accordance with the relationships (4) and (5) set forth above. 

When the inquiry of diamond 1750 is answered in the 
affirmative, diamond 1770 is entered, and inquiry is made as to 
whether the last level has been processed. If not, the block 
1730 is re-entered, to increment the level index, and the loop 
1780 is continued until processed values have been obtained for 
the nodes at all levels of the tree. 

Referring to Fig. 18, there is shown a flow diagram of the 
routine for computing the information cost of the nodes of the 
tree, so that an advantageous orthogonal basis can be selected. 
The block 1810 represents initializing the level index to the 
highest level (e.g., the last level in the illustration of Fig. 
3). The node index and the position index are initialized at 1 
(blocks 1815 and 1820). A node content count, which is used in 
the present embodiment to keep track of the number of processed 
values in a node that exceed a predetermined threshold, is 
initialized at zero, as represented by the block 1825. Inquiry 
is then made (diamond 1830) as to whether the value at the 
current position is less than a predetermined threshold value. 
If not, the node content count is incremented (block 1835), and 
the diamond 1840 is entered. If, however, the processed value at 
the current position is less than the threshold value, the 
diamond 1840 is entered directly. [At this point, the processed 
value could be set to zero prior to entry of diamond 1840 from 
the "yes" output branch of diamond 1830, but it is more efficient 
to handle this later.] Inquiry is then made (diamond 1840) as to 
whether the last position of the node has been reached. If not, 
the position index is incremented (block 1850), diamond 1830 is 
re-entered, and the loop 1855 is continued until all processed 
values of the current node have been considered. When this 
occurs, the node content count is stored for the current node (of 
the current level), as represented by the block 1860. Inquiry 
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is then made (diamond 1865) as to whether the last node of the 
level has been processed. If not, the block 1870 is entered, 
the node index is incremented, the block 1820 is re-entered, and 
the loop 1875 is continued until all nodes of the current level 
have been considered* Inquiry is then made (diamond 1880) as to 
whether the current level is the highest level. If so, there is 
no higher level against which comparison of parent-to-chi Idren 
node comparisons can be made. In such case, the level index is 
incremented (block 1885), block 1815 is re-entered, and the 
procedure just described is repeated to obtain and store node 
content counts for each node of the next-to-highest level. When 
this has been done, the inquiry of diamond 1880 will be answered 
in the negative, and the next routine (Fig. 19) will be operative 
to compare the level (which has just been processed to compute 
information cost of each node) with the children nodes of the 
previously processed higher level. 

In. particular, the level index is initialized (block 1905) 
to the highest level less one, and all nodes on the highest level 
are marked "kept". The node index is initialized (block 1910) 
and the node content count of the current node of the current 
level is compared (block 1920) to the sum of the node content 
counts of the two nodes which are children of the current node. 
[For example, if the current node is and the current level is 
L^, then the count for the current node is compared to the sum of 
the counts for nodes N^^^^^ and N^^ of level L^^^.] If the 
comparison shows that the parent has an equal or lower count, the 
parent is marked "kept", and the two children nodes are marked 
"not kept" (as represented by the block 1930). Conversely, if 
the comparison shows that the sum of two children nodes has a 
lower count than the parent node, each of the children nodes 
keeps its current mark, and the current parent node is marked 
"not kept" (as represented by the block 1940). In the case where 
the children nodes are preferred, the sum of the counts of the 
children nodes are attributed to the parent node (block 1945). 
By so doing, the lowest count will be used for subsequent 
comparisons as ancestors are examined. The attribution of the 
count to the parent node will not b problematic, since only 
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"kept" nodes will be considered in the next stage of processing. 
Inquiry is then made (diamond 1950) as to whether the last node 
of the current level has been reached. If not, the node index is 
incremented (block 1955)/ block 1920 is re-entered, and the loop 
1958 is continued until all nodes at the current level have been 
considered. Inquiry is then made (diamond 1960) as to whether 
the current level is level 1. If not, the level index is 
decremented (block 1965), block 1815 (Fig. 18) is re-entered, and 
the loop 1970 continues until all levels have been considered. 
At this point, the nodes which define the basis to be used have 
been marked "kept" [possibly together with some of their 
descendent nodes], and correspond, for example, to the shaded 
nodes of the Fig. 4 illustrations. 

Referring to Fig. 20, there is shown a flow diagram of a 
routine for generating output encoded words which, in the 
present embodiment, are collected in a frame which represents the 
encoded form of the data Xj^, Xg, X3 ... x^. For example, for an 
acoustical signal, the frame may represent a particular number of 
acoustical samples, for example 256 samples. As a further 
example, for a video signal, the freune may represent a frame of 
video, pori:ion thereof, or transformed frequency components 
thereof. The number of encoded words in a frame will generally 
vary as the information being encoded varies, and will also 
generally depend upon the level of the threshold employed, if a 
threshold is employed as in the present embodiment. Fig. 20 
illustrates an embodiment of a routine for generating a frame of 
words for the basis that was selected using the routines of Fig.s 
18 and 19 . A tree location index will be calculated which points 
to nodes in the tree in depth- first order (or so-called 
"pre-order") , as is well known in the art. The tree location 
'index is initialized to 1 at level 0, node 1 (block 2010). 
Inquiry is made (diamond 2015) as to whether the node at that 
tree location is market "kept", and, if not, diamond 2020 is 
entered directly, and inquiry is made as to whether the entire 
tree has been examined, as indicated by the tree location index. 
If the entire tree has been searched, block 2080 is entered and a 
"fram compl te" indication can be generated. If not, then loop 
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2011 is continued until a nod marked "kept" is encountered, or 
until the entire tree has been searched. If a node marked "kept" 
is encountered, block 2030 is entered, and the tree location 
index of this "kept" node is recorded in memory; suppose for 
example that it is called "X"* The position index in the node is 
initialized (block 2035). Inquiry is then made (diamond 2040) as 
to whether the value at the current position is above the 
predetermined threshold. if not, diamond 2055 is entered 
directly, and no word is generated for the value at the current 
position. If the value is above the threshold, block 2045 is 
entered, this block representing the generation of a word which 
includes the current level, node, and position, and the value at 
the position. The block 2050 is then entered, this block 
representing the loading of the just generated word into an 
output register. Inquiry is then made (diamond 2055) as to 
whether the last position in the current node has been reached. 
If not, the position is incremented (block 2060), diamond 2040 is 
re-entered, and the loop 2061 is continued until all positions in 
the node "X" have been considered. It will be understood that 
various formats can be used to represent the words. For 
example, a specific number of bits can be used for each of the 
level, node, position, and value. Alternatively, words could be 
of different length, e.g. depending on information content or 
entropy, with suitable delineation between words, as is known in 
the art. Also, if desired, all words in a particular node could 
be encoded with a single indication of level and node, with 
individual indications of position-value pairs. Inquiry is next 
made (diamond 2070), as to whether the last node location in the 
tree in depth-first order has been reached. If not, the tree 
location index is incremented (block 2070), and inquiry is made 
as to whether the new node is a descendant of "X", by a 
comparison of depth-first indices well known in the art. When 
this is the case, diamond 2065 is re-entered, and the loop 2071 
is continued until the first node which is not a descendant of 
"X" is encountered, or until there are no more nodes. When a 
first non-descendant of "X" is encountered, diamond 2015 is 
re-entered and the loop 2081 is continued until all nodes which 
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are both marked "kept" and have no ancestors mark d "kept" hav 
contributed to the frame. Such nodes contain a complete 
orthogonal group of wavelet-packet correlations (see also 
Appendix I, II/ and V). When either loop 2011 or the loop 2071 
terminates by exhaustion of the nodes, block 2080 is entered and 
a "frame complete" indication can be generated. If desired, the 
frame can then be read out of the encoder register. However, it 
will be understood that the encoder register can serve as a 
buffer from which words can be read-out synchronously or 
asynchronously, depending on the application. 

Referring to Fig. 21, there is shown a flow diagram of the 
decoder routine for processing frames of words and 
reconstructing the orthogonal basis indicated by the words of a 
frame* The block 2110 represents the reading in of the next 
frame. In the described embodiment, it is assumed that the 
frames of words are read into a buffer (e.g. associated with 
decoder processor subsystem 170 of Fig. 1), and the individual 
words processed sequentially by placement into appropriate 
addresses (which can be visualized as the selected basis nodes 
of a tree - as in Fig. 4), from which reconstruction is 
implemented in the manner to be described. However, it will be 
understood that individual words can be received synchronously or 
asynchronously, or could be output in parallel into respective 
tree nodes, if desired. Also, as was the case with the encoder, 
parallel processing or network processing could be employed to 
implement reconstruction, consistent with the principles hereof. 
In the routine of Fig. 21, the next word of the frame is read 
(block 2115), and a determination is made as to whether the node 
and level of the word is occurring for the first time (diamond 
2117). If so the node (and its level) is added to the list of 
nodes (block 2118). The value indicated in the word is stored 
(block 2120) at a memory location indicated by the level, node, 
and position specified in the word. It will be understood that 
memory need be allocated only for positions within the nodes 
designated by the read-in words. Inquiry is then made (diamond 
2130) as to whether the last word of the frame has been reached. 
If not, the block 2115 is re-entered, and the loop 2135 is 
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continued until all words of the frame have b en stored at 
appropriate locations. It will be understood that, if desired, 
the word locations (level, node, and position) could 
alternatively be stored, and the values subsequently recovered by 
pointing back to their original storage locations* 

During the next portion of the decoder routine, as shown in 
Fig. 22, the values in the nodes on the list are utilized to 
implement reconstruction as in the diagram of Fig.s 15A and 15B, 
with parent nodes being reconstructed from children nodes until 
the level zero information has been reconstructed. During this 
procedure, when a parent node is reconstructed from its children 
nodes, the parent node is added to the list of nodes, so that it 
will be used for subsequent reconstruction. This part of the 
procedure begins by initializing to the first node on the list 
(block 2210). Next, the block 2215 represents going to the 
memory location of the node and initializing to the first 
position in the node. inquiry is then made (diamond 2220) as to 
whether there is a non-zero value at the position* If not, 
diamond 2240 is entered directly. If so, the value at the 
position is mapped into the positions of the parent node, with 
accumulation, as described above in conjunction with 
relationships (10) and (11). Inquiry is then made (diamond 2240) 
as to whether the last position of the node has been reached. If 
not, the next position in the node is considered (block 2245), 
diamond 2220 is . re-entered, and the loop 2250 continues until all 
positions in the node have been considered. It will be 
understood that, if desired, a marker or vector can be used to 
indicate strings of blank positions in a node, or to point only 
to occupied positions, so that a relatively sparse node will be 
efficiently processed* In this regard, reference can be made to 
the abovereferenced U.S. Patent Application Serial No. 525,974. 
When the last position of the node has been considered, the node 
is removed from the list of nodes, as represented by block 2255, 
and inquiry is made (diamond 2260) as to whether the parent node 
is at level 0. If so, diamond 2270 is entered directly. If, 
however, the parent node is not at level 0, the parent node is 
added to the list of nodes (block 2265). Inquiry is then made 
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(diamond 2270) as to whether the last node on the list has been 
reached. If not, the next node on the list is considered (block 
2275), block 2215 is re-entered, and the loop 2280 is continued 
until processing is complete and the reconstructed values have 
been obtained. The decoder output values can then be read out 
(block 2290) . 

It will be understood that similar techniques can be 
employed at higher dimensions and in other forms (see e.g. 
Appendix V). More complicated tree structures, such as where a 
node has more than two children (e.g. Appendices II and V) can 
also be utilized. 

The invention has been described with reference to 
particular preferred embodiments, but variations within the 
spirit and scope of the invention will occur to those skilled in 
the art. For example, it will be recognized that the wavelet 
upon which the wavelet-packets are based can be changed as 
different parts of a signal are processed. Also, the samples can 
be processed as sliding windows instead of segments. 
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Appendix I 
Construction of Wavelet-Packets 

We introduce a new class of orthonormal bases of i^{R") by constructing a "library'" 
of modulated wave forms out of which various bases can be extracted. In particulai\ 
the wavelet basis, the w^h functions, and rapidly oscillating wave packet bases axe 
obtained. 

We'll use the notation and terminology of [D], whose results we shall assimie. 

§1. We are given an exact quadrature mirror filter h{n) satisfying the conditions of 
Theorem (3.6) in [D], p, 964,-1.6. 

53 h{n - 2k)h{7i - 2£) = Sk^e , ^ h{n) = 

We let gk = hk+i{-l)^ and define the operations Fi on i^{Z) into ^£'^{2Z)'' 

Fi {sk}(i) = 2 53 Sk9k^2i- 

The map F{sk) = Fo{Sk) © Fi(sfc) € £^(2Z) © i^{2Z) is orthogonal and 

(1.1) F^Fo + F^Fx=I 

We now define the following sequence of functions. 

(12) / = ^^5: hkWr,{2x - k) 

I W^2n+i (x) = E gkW„(2x -k)' 

Clearly the function Wb(x) can be identified with the function ip in [D] and Wi with 
the function ip. 

Let us define moiO = 75 S Afce"'*^ and 
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Remark: The quadrature mirror condition on the operation F = (i^o? ) is equiv- 
alent to the unitarity of the matrix 


I^Jdng Fourier transform of (1.2).when n = 0 we get 
Wb(0 = mo(f/2)Wb(^/2) 


I.e., 

oo 

and 

WiiO = mi(e/2)IVb{e/2) = miU/2)moUM)moU/2^) ■ • • 
More generally, the relations (1.2) are equivalent to 

(1-3) w^{o = ^rn,jmn 

and n = 2 €j2^-^ {sj = 0 or 1). 
We can rewrite (1.1) as follows. 

(1.4) W^2„(a: - ^) = %^ J3 hj.2tW„(2x - j) = Fo{W„(2x - j)}(e) 

W2n+i{x -£) = V2j29i-2eWr,{2x-j) = J'i{W^„(2x - j)}(£) 

• where Wni2x — j) is viewed as a sequence in j for (s,n) fixed. Using (1.1) we find: 
(1.5) 


W„ix - j) = y/2^hj.2iW2^ (f - *) 9j-2iW2,^i (I - i) 


In the case n = 0 we obtain: 
(1.6) Woix -k) = y/2j2 hk-2iWo (| - ») + gt-aWr (| - ») 
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from which we deduce the usual decomposition of a function / in the space SIq (Vq in 
[D]) i.e., a function / of the form 

More generally, if we define 

(1-7) nn = {f:f = J2uJtW„{x-k)}. 

We find 

(1.8) /(X) = J2 ^Foio^k){i)W2„ (f - 0 + E ^FiMii)W2.+i (f - i) 
or 

y/2f{2x) ^ h + g € fi2n 5 € ft2».+i 

We now prove 

Theorem (1.1). The functions Wn{x — k) form an orthonormal basis of L^{Ii), 

Proof, We proceed by induction on n, assuming that Wn{x — k) form an orthonormal 
set of functions and, proving that, W2rt{x — k), W2n^i{x - fc) form an o.n set. 

By assumption || V2/(2a;)||| = 53 a/| if / € ^U from the quadrature mirror condition 
on {Fo^Ft) we get 

Since Fo{<ifk){i) = M», Fo(a;jb)(i) = i/^ can be chosen as two arbitrary sequences of i'^ 
(arising from a; = F^fii + F^t/i) it follows that 

/ 1 E A^W^2n(x - 0 + X) ''*W^2.+i - OP = 13 + 53 X/? 

which is equivalent to W2n{x - Oi W2n+i(a: - j) being an o.n set of functions. 

Let us now define Sf = ^/2/(2x). Formula (1.8) shows that = fign 0 ^211+1 as 
. an orthogonal stmi or, 

(1.9) SQq - tto = J^i 

5^no - S^Slo = *f22 ® = ^4 ® fis ® fie e ^7 or 
^^'Jio - ^^'"^lio = figt-t ® n2fc-i+i ' • ' ® Sl2t'-l 

and 
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s^Oo = J2o ® Hi e • • • e ^2^-1 

More generally, we let Wk = S'^'^^SIq - S^SIq = ^^^i = S'^Wi. Therefore we have 
Proposition (1.1). 

>Vib = 5*^1^1 = f22* e ^2^+1 © • • • © J^2fc+i-l- 

Alternatively, the functions 

W^ix^j) i G Z 2*^ < n < 2^-+^ 

form an ortbonormal basis ofWk. 

Since the spaces Wk are mutuaUy orthogonal and span Ir^{R) see [D], it follows that 
Wn{x — j) are complete. 

§2. Orthonormal bases extracted out of the "library" 2^^^Wr^{2^x - k). 

We start by observing that 2^^^Wii2^x - A:) form a basis of Wj as we vary A: and a 
basis of L^{It) as j, k vary. This is the wavelet basis constructed in p]. The following 
simple generalization is useful to obtain a better localization in frequency space. 

Theorem (2.1). The functions 

2^^^Wni2^x - fc) 7 = 0, ±1, . . . , fe = 0, ±1, ±2 

2^ < n < 2^+^ for Sxed £ form an ortbonormal basis of L^{R) ^ = 0, 1, 2, 

Remark: This is a wavelet basis with dyadic dilations and 2^ fundamental wavelets. 

Proof. We have seen, in Proposition (1.1), that Wn{x - fc) 2^ < n < 2^+^/c € Z form 
an o.n basis of We, from which we deduce that 2^/^Wrt{2^x - k) form an o.n. basis of 
Wi^j spannmg L^(R) for j € Z fe € Z. 

In the next example we vary j and n simultaneously to obtain a basis whose number 
of oscillation is inversely proportional to the length of its support. 

Theorem (2.2). Let t{n) = [loggn] i.e., 2«(") <n< 2^+^ tben 

TV'„(2^(">x - A:)2^^">/2, Wn{2^^"^+'x - k)2^^ 
form an o.n basis of i^(R). 

Pn>o/. Fix £(n) = ^, and consider n with 2^ < n < 2^+^ as seen in the proof of 
Theorem (1.2) 

Wn{2^3: - fc)2^/^ form an o.n basis of W2e 
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and 

Wn(2^+^x-A:)2^^+i)/2 form an o.n basis of W2£+i 

since = ^ e>V2£ ® J2 ®W2£+i, we have a complete basis. 
These can be generalized as follows. 

Theorem (2,3) • Let sl collection be given such that the dyadic intervals Ie,n = 

[2V 2^(n+l)) form a disjoint covenng of (0, oo), then 2^/^Wn{2^x-k) form a complete 
ortbonormal basis of L^(R), 

This theorem becomes obvious in the following case. Let 

lO f<|e|<7r 
a periodic function of period 2ir, and mi(^) = 1 - mo(^). Let 

3=1 

then 

»»<|e/'r| <n + l 
1 0 elsewhere 
and the orthonormal basis ufn{x — fe) in Fourier space is 

e*ed>n(e) 

wMch is the simplest variation on a "Tvindowed" (2 windows) Fourier transform. The- 
orem (1.4) is obvious in this case. This theorem is also easy to understand from the 
point of view of subband coding as we shall see. 

§3. Subband coding and expansions in terms of Wn 

We assume, given a function which, on the scale 2~*^, is well approximated as 

(3.1) /(X) = 53 slWo{2^x - A:)2^/2 

as seen in (1.8) 
(3.2) 
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The coefficients of /o £tre given by Fois^)- 
The coefficients of fx are given by i^i(s°). 
Continuing an application of (1.8) gives 

= ^ (f 2") + (f 2") + A. (f 2") + /„ (^22") } 

where /oo>/io axe obtained by decomposing /o 

/019/11 2^ obtained by decomposing /i 

foo € flo>/io € f2i,/oi € ft2i/ii € fla- If we continue this decomposition and observe 

that the binary tree corresponds to the realization of n as n = ^€j2^~^ and that after 

1 

£ iterations we get 

(3.3) fix) = 2^^-^>/2 J2 /«(^2^-^) with /„ € 
and U(x) = W^T.(2^-^x - k))Wr,(2^''^x - fc)2^"^ with 

k 

(3.4) 2:^(/, PV„(2^-« - k)) = f-^.Fe, . . . JV, {^2} 

We therefore obtain a fast 2^N algorithm to calculate all coefficients for "all functions 
in our library*' for scales —N < j < 0, The procedure is analogous to subband coding. 

§4. EUgher dimensional libraries and bases 

We define the higher dimensional wavelet basis as follows: Let 

Wi = spsin{Wo{xi-kj)Wi{x2'k2hWi(xi^kx)WQ{x2-k2lWi(^^ 

Wk = S'^Wi where Sf (X1X2) = 2/(2xi,2x2). Clearly 

The basic 2-dimensional ^'Ubrary" consists of all functions obtained as tensor prod- 
ucts of the one dimensional library i.e., 

2iO*i+ia)|^^^(2iia;i - ki)Wn^{2^^X2 - 62). 

We will restrict our attention to the sublibraries obtained by dilating both variables 
by the same dilations although the case j2 = rjj r fixed is of independent interest. 
The two dimensional basis corresponding to Theorem (2.1) is given in 
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Theorem (4,1). Fix i, then for (fci, fca) € j e Z 2^ < n < 

2^^^W^{2^xj - fei)2^Wo(2^+^X2 - fca) 

2^Wo{2^'^^xi - fci)2J/2w^(2Jx2 - &2) 

2^ Wn, - fci ) Wn, (2^X2 - fc2 ) 2^ < ni < 2^+^ 

form an ortlzonormaZ basis of . 

Theorem (4.2). Let (n> = <ni,n2> ^ 2™«(^("»)^<"'» where 
i{ni) = [log2 ni) ni > 0 n2 > 0 tien for (fci, ^2) € Z^, 

<2n>W„, (2<n)a:x - fci )W^^(2(n)x2 - fe2) 

form an ortltonormaZ basis (wavelet packet basis of L^(R?)). 
Proof. Assume £ = max (^(ni)£(n2)) = ^(^2), i.e. ni < 7x2. Then 

2^Wr,, (2^X1 - fei ) Wn, (2^X2 - ^2) 

for 

0 < ni < 2^ 2^ < 712 < 2^+^ 

spaji (using 1.9) Proposition (1.1) S^^SIq^^ 0 5^ni,x2 i'®- subspace spanned by 
Wo(2^^Xi — fci)Wi(2^^X2 — fc2)* Ck)nsideration of the other 2 cases yields a total span 
of W2£, similarly we obtain >V2£+i and the theorem is proved. 

Reference 

[D] Ingrid Daubechies, Orihonormcd bctsea of compactly supported wavelets^ CommuDications on 
Pure and Applied Mathematics XI«I 909-966. 
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Appendix II 
Best- Adapted Wavelet-Packet Bases 

Introduction. By using ertra filters, it is possible to introduce fast wave packet 
transformations which decimate by arbitrary numbers. Such transformations generalize 
algorithms which decimate by 2. The method produces new Ubrari^ of orthonormal 
basis vectors. We introduce an algorithm for selecting a most efficient representation 
firom the library^ and prove that its complexity is 0{N) for a sequence of length N. 
We discuss some of the analytic properties and applications of such representations. 

Aperiodic filters and bases in Consider first the construction of bases on P. 
Let p be a positive integer and introduce p absolutely summable sequences /o, . . . , /p—i 
satisfying the properties: 

(1) For some £ > 0, X)^ \fi{m)\\fn\^ < oo, 

(2) Em/o(^) = l»fori = 0,l,...,p-l,and 

(3) fi{''^)fj{''^ + kp) — Si^jSk, where 6 is the Kronecker symbol. 

To these sequences are associated p convolution operators J^o^ • • • ^ -^p— i aJ"! their ad- 
joints i^Q , . . . , -Fp_i defined by 

Fi:l^-^l\ FMk) = 53 fi{7n+pk)v{mh 

tn 

Ftie^i", F:v{m)=^YlMm^pkMk). 

k 

These convolution operators will be called Slters by analogy with quadrature mirror 
filters in the case p = 2. They have the foUowing properties: 

• Lemma. For i^j = 0, 1, . . . , p — 1, 

(1) FiF^^O, ifi^j. 

(2) FiFt = /, 

(3) F*Fi is an orthogonal projection of , and for i ^ j tlie ranges of F*Fi and 
FJFj are orthogonal, and 

(4) Fo*i^o+ •* + F;LiFp-i =/. 

Research supported in part by ONR Grant N00014-88-K0020 
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Proof. Properties (1) and (2) follow by interchanging the order of integration: 
FiPrvik') = 53 X) ^i^"^ + pk')Mm+pk)v{k) 

m k 

= S ( E M^')Mm' ^p[k - k'])) v{k) 
k \rn' / 

= E Si^ySk^k' v{k) 
k 

^ft;(&'), if*=i, 
lO, ifi#j, 

by the orthogonality properties of /». 

Property (3) follows from (1) and (2): F^FiF^Fi = F^Fi, and F^FiF^Fy = 0. 
Orthogonality is easily shown by transposition. 

To prove (4), let mj(f) = T^kfji^)^^ ^® (bounded, Holder continuous, peri- 
odic) function detemuned by the filter for j = 0, . . . ,p - 1. Then fy{k) = mj{k) 
IS a real number, and each F^Fi is unitarily equivalent to multiplication by \mi\^ on 
£2(«7r,7r). 

Now PlancherePs theorem gives 



In particular, \mj\^ has integral 1, and the Fourier coefficient (|mjp)r(Zp) vanishes if 
1^0. This is equivalent to the average of Imy over ^-f27r/p, . . . , ^+27r(p— 1) /p} 
being identically 1. 

The same vanishing is true of the Fourier coefficients of the cross terms my my, and 
for those it also holds when 1=0. Thus, the average of rnj{^)mjr(^) over + 
27r/p, . . . ,f -h 2^{p — l)/p} vanishes identically. Hence, the conditions on the filters fi 
are equivalent to the unitarity of the matrix 

/ mo(0 mo(e + ^) ... mo(e + 2£le=ii) \ 

But then Sj^J |mfc(e)|2 = 1 for aU Thus F^Fo + ■■■+ F;_^Fj^i is unitarily 
equivalent to multiplication by 1 in L^(—ir,ir), proving (4). □ 

With this proposition we can decompose into mutually orthogonal subspaces 
Wo* ± ••■ ± W^i, where Wl = F;Fi{l^) for i = 0 The map F. finds 
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the coordinates of a vector with respect to an orthonormal basis of W/, Since each 
FiW^ = Fi{P) is another copy of there is nothing to prevent us from reapplying the 
filter convolutions recursively. At the mth stage, we obtain = ± - - ± PF^„i, 
where W;^ = F*^ . . . F^i^Fn^ ...Fn^l^) and . . . ni is the radix-j? representation of 
n. The map Fn„ . . . Fy^ transforms into standard coordinates in W^. For convenience, 
we will introduce the notations F^.= . . . F^, and F^* = F*^ . .^F^, 

The subspaces WJP^ form ap-ary tree. Every node is a parent with p daughters 
• • * > The root of the tree is the original space which we may label 

W§ for consistency. CaU the whole tree W. 

Now fix m and suppose w belongs to TT^, where 0 < n < p"* - 1, and F^w = Ck is 
the elementary sequence with 1 in the fcth position and O's elsewhere. The collection 
of all such w forms an orthonormal basis of with some remarkable properties. In 
particular, if p = 2 and the filters Fq and Fi are taken as low-pass and high-pass 

quadrature mirror filters, respectively, then the spaces Wg"' W^^^ are all the 

subbands at level m. These have been used for a long time in digital signal processing 
and compression. In an earher paper we described experiments with an algorithm for 
choosing m so as to reduce the bit rate of digitized acoustic signal transmission. This 
produced good signal quality at rather low bit rates. 

The tree contains other orthogonal bases of Wg. In fact, it forms a library of bases 
which may be adapted to classes of functions. The tree structure allows the library to 
be searched eflBciently for the extremum of any additive functional. 

To every node in W we associate the subtree oT all its descendants. Define a graph 
to be any subset of the nodes of W with the property that the union of the associated 
subtrees is disjoint and contains a complete level Wg", . . . , for some m. The 

singleton {W§} is a graph, for example, with m = 0. The following may be called the 
graph theorem. 

Theorem. Every graph corresponds to a decomposition of l^ into a 6nite direct sum. 

Proof. Every graph is a finite set, of cardinality no more than p"* for the m in the 
definition. Fix a graph, and suppose that W^^^ and WJ^^ are subspaces corresponding 
to two nodes. Without loss, suppose that mj < m2. Then W;^^ is cont^ed in a 
subspace for some n ni. Since the subspaces at a given level are orthogonal, 

we conclude that W^^ ± W^^^' . 

To show that the decomposition is complete, observe that a node contains the sum 
of its daughters. By induction, it contains the sum of all of the nodes in its subtree. 
Hence a graph c ntains the sum of all the subspace at some level m. But this sum is 
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allofZ2. □ 

Theorem. Gr&pbs are in one-to-one correspondence with Snite disjoint covers of [0, 1) 
byp-adic intervals Ii;^ = p""*[n,n + 1), n = 0, 1, . . - ,p"* - 1. 

Proof. The correspondence is evidently WJP^ iJJ*. The subtree associated to 
corresponds to all p-adic subintervals of I^, The details are left to the reader. □ 

Analytic properties of graphs: continuous wave packets. Each filter Fy- (and 
its adjoint Fj) maps the class of exponentially decreasing sequences to itself. Likewise, 
the projections FJJ**F^ preserve that class. In practice, we shall consider only finite 
sequences in P. For actual computations the filters must be finitely supported as well. 
Convolution with such filters preserves the property of finite support. Let the support 
width of the filters be r, and let be the maximum width of any vector of the form 
jPj* . . . F^^ (cfc). Then zq = 1 and Zm+i = pZm + r — p. By induction, we see that 

Coifhian and Meyer [CM] have observed that the basis elements Fl^*eis are related 
to wave packets over R. A slightly generahzed paraphrase of their construction follows. 
Many of the basic facts used were proved by Daubechies in [D]. 

Let w be a function defined by w{^) = Hy^i "^(f/p' )' where mo is the analytic 
function defined by Fq^ ss above. Then w has mass 1, decreases rapidly, and is Holder 
continuous, as proved in [D]. If mo is a trigonometric polynomial of degree then 
w is supported in the interval [— r,r]. Arranging that w has r continuous derivatives 
requires mo with degree at most 0(r). See [D] for a discussion of the constant in this 
relation for p = 2. Put = to, and define the family of wave packets recursively by 
the formula w^^j{t) = J^-oo "~ «)- This produces one function for each 

pair (m, n), where m = 0, 1, . , . and n = 0, 1, . . . ,p"* — 1. 

We can renormalize the wave packets to a fixed scale p^. Write 

Then Wq q^^ is a collection of orthonormal functions of mass 1, concentrated in intervals 
of size 0{p~^), This makes them suitable for sampling continuous functions. Let x(t) 
be any continuous function, and put 

J — oc 

We may use Sq(&) as a representative value of x{t) in the interval 7^ = p~^[fc, k + 1). 
The closeness of the approximation to valu^ of x depends, of course, on the smoothness 
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of X. Suppose that x is Holder continuous with exponent e. Then if to is any point in 
/jf*, we have 

|x(to) - sl{k)\ = I / {x{to) - x(t)) tog(p^t - k)dt\ = 0(p-^^). 
Jit 

We can also take advantage of differentiability of x if we construct with vanish- 
ing moments. Given d vanishing moments and d derivatives of x, the approximation 
improves to |x(to) - s^{k)\ = 0(p-''^). 

The map x h-^ fig sends L^(R) to and pulls back the orthonormal bases of con- 
structed in the last section. To see this, define 8^{k) = (x^w^^^^j,). By interchanging 
the order of recurrence relation and inner product, we obtain the formula = FS^sg. 
Thus, the coordinates s^ik) are coefficients with respect to an orthonormal basis of 

The resulting subspaces of ^^(R) form a finer type of multiresolution decomposition 
than that of M^at [Ma]. The coordinates s'^(k) are rapidly computable. They contain 
a mixture of location and frequency information about x. 

Ordering the basis elements. The parameters n,m,fc,X in w^ „,j, have a natural 
interpretation as frequency, scale, position, and resolution, respectively. However, n is 
not monotonic with frequency, because our construction yields wave paxdcets in the so- 
caQed PaJey (natural, or p-adic) ordering. The following results show how to permute 
n H-» n' into a frequency-based ordering. 

Theorem. We can choose rapidJy decreasing filters Fq - 1 such that u^n.m.fc 

is concentrated near the interval /f'""*, and w^:,^,^ is concentrated near the interval 
J^, where n n' is a permutation of the integers. 

Proof. The first part is evident. For any family of exponentiaUy decreasing filters, 
tyg decreases exponentially away from [0,1). ti;^,„,,fc is its dilate and translate to the 
interval J^^"*. Likewise, w^„,j, has the same concentration as ti;^,^ ^' ^ 
filters Ft are uniformly exponentially decreasing. 

The second part follows from the Fourier transform of the recurrence relation: 

= (p'''^Mk)e-'^^^^^ tS-U/p) = p-'m,(e/p)t&-(e/p), 

where mj is the multiplier defined above. Recall that Xly^o = 1 and that 

Tno(O) = 1. Thus, the periodic functions |mjp form a pai'tition of unity into p functions. 

with 0 being in the support of mo alone. 
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Now suppose for simplicity that we have chosen filters in such a way that |t71j(^)| = 
]Cte-oo Xi^Ii j+i) (f — 2nk). Such rrij may be approximated in (— tt, tt) as closely as 
we like by multipliers arising from exponentially decreasing filters. In this simple case, 
it is inunediate that Wq(() = wio(f/p)I(-ir,x) is the characteristic function of (— 7r,7r), 
so that t&o'.o.o is characteristic function of (— 7rp'^,7rp^), Likewise, t&^i,o is the 
characteristic function of — 1, — j] U "rrp^^ljyj + 1). 5Vom the recurrence 

relation, we see that Wn,m,o characteristic function of the miion of the 

intervals ±7rp^"*[n',n' + 1), where n k+ n' is a permutation. These intervals cover 
p^(— 7r,.7r) as n = 0, . . . — 1. The permutation n n' is given by the recurrence 
relation 


•r n n / _l m/ / ^P-^^^ if n' IS 

if n = 0,...,p-l; {np+jY=< , . . 

L n'p + (p — 1) — J, if n' IS 


even. 


odd. 

Write fij for the jth digit of n in radix p, numbering from the least significant. Set 
Um = 0 if n has fewer than m digits. Then the recurrence relation implies that 
Uj — 7r(n^-^i,n^), where 


, . f if X is even, 

Ip— 1— y, ifxis odd.. 


For each value of the first variable, tt is a permutation of the set {0, . . . ,p — 1} in the 
second variable. Thus the map n' i— ► n and its inverse ni-* n' are permutations of the 
integers. It is not hard to see that these are permutations of order 2 if p happens to be 
odd. Otherwise they have infinite order, as may be seen by considering an increasing 
sequence of integers n' all of which have only odd digits in radix p. □ 

Corollary. With Slters Fq, . . . ,i^p-i chosen as above, we can modify the recurrence 
relation for toj^'^m^ such that w^^^j. is concentrated near the interval I^. 

Proof. Simply reorder the functions by using the alternative recurrence relation: 


m+i ^ / Sfc fj(k)w^{pt - A:), if n is even, 

V+i^ ) \ /p_l.j(A:)«;-(pf - fc), if n is odd. 


Since we are enforcing n = n' at each level m, we are composing with the permutation 
defined above. Of course, this algorithm has complexity identical to the original. □ 

Periodic filters and bases for R''. A sampled periodic function may be represented 
as a vector in R** for some d. In this case let p be any factor of d. Introduce as filters 
a family of p vectors {fi 6 R** :i = 0 p — 1}. These are obviously summable. 
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Suppose in addition that they are orthogonal as periodic discrete functions, i.e., that 

E^i/iM/i("» + *=P mod d) = Si-jSk. 

Let the associated convolution operators be {Fq, . . . , Fp-i}, defined as above by 


p.^ . _» Rd/p, Fiv{k) = 5Z d)v{m), for fc = 1, 2, . . . , d/p, 

m=l • 
d/p 

: TL'^" -* B*, F*w{m) = ^ fi{m + pk mod d)v(k), for m = 1, 2, . . . , d. 

fc=i 

The reduction modulo d is IntentionaUy emphasized. These operators satisfy conditions 
similiar to those of aperiodic filters: 

Proposition. 

(1) FiF; = 0, if* ^3, 

(2) FiFt = Ii/p 

(3) F^Fiisa. rank d/p orthogonal projection on R**, and for t # j the ranges of 
FiFi and F^Fj are orthogonal, 

(4) FSFo + • • • + F;_jFp_i = Id 
where Id is the identity on R'*. 

Proof. The proof is nearly identical with the one in the aperiodic case. □ 

The decomposition suggested by equation (4) may be recursively applied to the p 
subspaces R-*/' by using additional filter families. For d = pj . . .pi and 0 < n < d. 
we have a unique representation n = m + naPi + ngPaPi + • • • + ntPL-i • • • Pi , where 

0<ni<Pi- This defines a one-to-one correspondence between {0 d - 1} and an 

index set of L-tuples J = { (m , .... n^) : 0 < < Pi}. We caai construct a b^is of R" 

whose elements are indexed by /. Forn = (nj fij,) € /, define = F^^ . . . F„, , 

where F» is a famUy of Pi periodic filters. Then F^:*F^ is an orthogonal projection 
onto a 1-dimensional subspace of R''. This is shown by induction on the rank in (3). 
'Now let Wj^ be the range of this projection. The coUection {u„ = F^*l : n € J} 
of standard" basis vectors of W;^ will be an orthonormal basis of R", and the map 
pi . R_d _» R gives the component in the tin direction. 

"as before, we are not limited to the basis defined by the index set /. Products of 
fewer than L filters form orthogonal projections onto a tree of subspax:es of R". A node 
arising from a product of m filters wiU correspond to the subspace WT = F™*F;rR'*. 
where n = n, -f- ■ • • + n^Pm_i • • - Pi indexes a composition of m filters. The tree wiU 
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be nonhomogeneous in general, although all nodes i levels from the root will have the 
same number pi of daughters. Define a nonbomog^eneous graph as a finite union of 
nodes whose associated subtrees form a disjoint cover of some level m < L, A graph 
theorem holds for this tree of subspaces as well. It and its corollary may be stated as 
follows: 

Theorem. Every nonbomogeneous graph corresponds to an orthogonal decomposition 

Corollary. Graphs are in one-to-one correspondence with Snite disjoint covers of [0, 1 ) 
by intervals of the form /JJ* = (pi . . .pm)~^[n, n + 1). □ 

Any permutation of the prime factors of d gives a (possibly different) basis. 

Smooth, filters. Some filter sequences have a smoothness property: 

Definition. A summable sequence f is a smooth Slter (of rankp) if there is a nonzero 
solution 4> in i^(R) fl i^(R) n C~{R) to the functional equation 

<f>{x) = 53 f{m)4>{px + m). 

m 

A filter will be said to have smoothness degree r if it satisfies this definition with 
replaced by C. Daubechi^ has shown in [D] that finitely supported filters of any 
degree of smoothness may be constructed in the case p = 2. An obvious consequence 
is that smooth filters exist in the case p = 2'. For arbitrary p, we may construct a 
filter fanuly as above subject to additional constraints. 

A continuous solution to the functional equation (3) always exists for a sequence 
/ satisfying the three conditions at the top of this article. Its Fourier transform may 
be constructed by iteration: 

oc 

ka = ko)l[rn{ap''h 

where m(^) = p"^/^ /(A:)e~*^^ is the multipUer con*esponding to the filter convo- 
lution in the integral equation. If ^ is nonzero, then ^(0) 7^ 0. so it may be assumed 
that ^(0) = 1. Now the sequence {/(A;)|fc|*} converges absolutely, so m(^) is Holder 
continuous of degree e. But also, m(0) = p^/^ ~ ^' ^ 0 the 

estimate |m(^) — 1| < C|f |^ holds. This implies that the infinite product converges. 
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Appendix III 
Nonstandard Matrix Multiplication 

Wave packets. Define wave packets over in the usual way. For a pair P = {p* } , Q = 
{qi) of quadrature mirror filters (QMFs) satisfying the orthogonality and decay con- 
ditions stated in [CW], there is a unique solution to the functional equation 

= V^53pj^(2t-i). 

i€Z 

Put w = itio,o,o = <f>^ and define recursively 

i€Z 

ti'2n+l,0.0(t) = 9i^n.0,0(2t - j). 

Then set t«„mfc(t) = 2^I^Wr^m{2'^t - fc). Write W(R) = {to^mfc : n,m,fc € Z} for the 
collection of functions so defined, which we shall call wave packets. 

The quadrature mirror filters Q may be chosen so that >V(R) is dense in many 
common function spaces. "VWth the minimal hypotheses of [CW], >V(R) will be dense 
in i2(R). Using the Haar filters P = {l/>/2, 1/V2}, Q = {l/v/2, -1/a/2} produces 
W{R) which is dense in jCfP(R) for 1 < p < oo. Longer filters can generate smoother 
wave packets, so we can also produce dense subsets of Sobolev spaces^ etc. 

Basis subsets. Define a basis subset a of the set of indices {(n,m, fc) € Z^} to be 
any subcollection with the property that {tOnmfc : (n,m,fc) € a} is a Hilbert basis for 
2/^(R). We characterize basis subsets in [Wl]. Abusing notation, we shall also refer to 
the collection of wave packets {tWnmt : (n,m,&) 6 a} as a basis subset. 

Since L*^ rxU" is dense in Zi^ for 1 < p < oo, with certain QMFs a basis subset will 
also be a basis for U*. Likewise, for nice enough QMFs, it will be a Hilbert basis for 
the various Sobolev spaces. 

Since L^(R)®I/^(R) is dense in L^(R^), the collection of vectors {wx®wy : v)x ^ 
W{X)^WY ^ is dense in the space of Hilbert-Sdiniidt operators. Call a C a 
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basis subset if the coHection {Wnxmxkx ®«'nv»nyfcr '• fcjf,ny,n»r, fri') e a} 

forms a Hilbert basis. Such two-dimensional basis subsets are characterized in [W2]. 

Ordering wave packets. Wave packets Wnmk can be totally ordered. We say that 
wKw'U (m,n, k) < (m',n', fc')- The triplets are compared lexicographicaUy. counting 
the scale parameter m as most significant. 

Tensor products of wave packets inherit this total order. Write wx = Wn vmxfc.v 
etc. Then we wiU say that wx < w'x^w'y wx < w'x, or else if twx = w'x but 

wy < vj'y. This is equivalent to 

(m.x,nx,kx,mY,nY,W) < {mx,n'x,k'x,mY,nY,kY), 

comparing lexicographically from left to right. 

Define the adjoint order <* by exchanging X and 1" indices, namely wx ® t^v <* 
w'x^v/'y if and only if toy 0 wx <* t^y ® -^'x- This is also a total order. 
Projections. Let denote the space of bounded sequences indexed by the three 
wave packet indices n, m, k. With the ordering above, we obtain a natural isomorphism 
between l°° and W^. There is also a natural injection : Z-^CR) '-* given by 
Cr^mk = (w^tfnm*), for V € i*(R) and w^rru. € W(R). If a is a basis subset, then 
the composition J* of J* with projection onto the subsequences indexed by <t is also 
mjective. is an isomorphism of ^^(R) onto l^{<r), which is defined to be the square 
summable sequences of whose indices belong to a. 

We also have a map R} -* L^(R-) defined by 

(«.Tn.,fc)6ZS 

This map is defined and bounded on the closed subspace of isomorphic to under 
the natural isomorphism mentioned above. In particular, is defined and bounded 
on the range of for every basis subset a. The related restriction Rl : -* 
LHR) defined by Ric(t) = E(n,m.fc)€a CnmkW„^{t) is a left inverse for and In 
addition, J^Rl is a projection of Likewise, if XliO^ = 1 and Rl^ is one of the 
above maps for each i, then Ei «i-Ra, is also a projection of WMt is an orthogonal 
projection on any finite subset of W^- 

SimUarly, writing for x W^ the ordering of tensor products gives a natu- 
ral isomorphism between Z~ and W^. The space L^B?). i.e.. the Hilbert-Schmidt 
operators, inject into this sequence space in the obvious way, namely M ^ 


SUBSTITUTE SHEET 


wo 91/18361 


PCr/US91/03504 


- 39 - 

in-3 

{M^Wnxmxhx ® ^nymj-fcy)- Call tUs injection J^. If <7 is a basis subset of W^. 
then the composition of with projection onto subsequences indexed by a is also 
injective. is an isomorphism of L^{R?) onto /^(a), the square summable sequences 
of whose indices belong to a. 

The map : — > L^(R^) given by R?c{x^y) = cxi^ti;x(x)t(;v(y). is bounded 
on that subset of naturally isomorphic to i^. In particular, it is bounded on the 
range of for every basis subset a. 

We may also define the restrictions of R^ to subsequences indexed by a, defined 
by Rlc{x,y) = J2 {-ujx *vn-)^<T ^^y'^^^^^'^^^y^' There is one for each basis subset a of 
W^, Then ilj is a left inverse of and J^, and J^R^. is a projection of W^. As before, 
if CKt = 1 and is a basis subset for each i, then R^^ is also a projection of 

>V^. It is an orthogonal projection on any finite subset of W^, 

Applying operators to vectors. For definiteness, let X and Y be two named copies 
of R. Let V € L^{X) be a vector, whose coordinates with respect to wave packets form 
the sequence J^v = {(v^wx) • tax € W{X)}, 

Let M : L^{X) — ♦ L^{Y) be a Hilbert-Schmidt operator. Its matrix coefficients 
with respect to the complete set of tensor products of wave packets form the sequence 
J^M = {{M,wx <S> wy) : wx e W{X),wy e >V(i^)}. We obtain the identity 

WA-€W(X) 

This identity generalizes to a linear action of on defined by 

C{v)nmk = Cnmkn'TTt'k'Vn'm'k'- 
(n'm'Jb') 

Now, images of operators form a proper submanifold of W^. Likewise, images of vectors 
form a submanifold W^. We can lift the action of M on v to these larger spaces via 
the commutative diagram 

-I 

The significance of this lift is that by a suitable choice of a we can reduce the complexity 
of the map J^M, and therefore the complexity of the operator appUcatiou. 


M 
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Composing operators. Let X, Z be three named copies of R. Suppose that 
M : L^{X) L^(Y) and N : Ir^(F) L^{Z) are Hilbert-Schmidt operators. We 
have the identity 

This generalizes to an action of on W^, which is defined by the formula 

n"m"Jb" 

where c and d are sequences in W^, Using J^, we can lift multiplication by TV to an 
action on these larger spaces via the commutative diagram 

LHB?) > L^iB?) 

N 

Again, by a suitable choice of a the complexity of the operation may be reduced to 
below that of ordinary operator composition. 

Operation counts: transforming a vector. Suppose that M is a non-sparse opera- 
tor of rank r. Ordinary multiplication of a vector by M takes at least O(r^) operations, 
with the minimum achievable only by representing M as a matrix with respect to the 
bases of its r-dimensional domain and range. 

On the other hand, the injection will require ©(r^pogrj^) operations, and each 
of and B} require O(rlogr) operations. For a fixed basis subset a of W^, the 
appUcation of J^M to J^v requires at most #|JjM| operations, where #|?7| denotes 
the number of nonzero coefficients in 17. We may choose our wavelet library so that 
#|J^Af| = O(r^). Thus the multiplication method described above costs an initial 
' investment of 0(r^[logr]^), plus at most an additional 0{r^) per right-hand side. Thus 
the method has asymptotic complexity 0{r^) per vector in its exact form, as expected 
for what is essentially multiplication by a conjugated matrix. 

We can obtain lower complexity if we take into account the finite accuracy of our 
calculation. Given a fixed matrix of coefficients C, write Cs for the same matrix with 
all coefficients set to 0 whose absolute values are less than 6. By the contmuity of the 
Hilbert-Schmidt norm, for every e > 0 there is a 5 > 0 sucli that \\C-Cs\\hs < Given 
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M and e as well as a library of wave packets, we can choose a basis subset a C 
so as to minimize H=\{J^M)s\, The choice algorithm has complexity 0(r^[logr]^), as 
shown in [W2]. For a certain class of operators, there is a library of wave packets sucli 
that for every fixed 5 > 0 we have 

(S) #|(J^M)fi| = 0(rlogr), 

with the constant depending, of course, on 6, We will characterize this class, give 
examples of members, and find useful sufficient conditions for membership in it. For the 
moment, call this class with property S the sparsiBsLble Hilbert-Schmidt operators S. 
By the estimate above, finite-precision multipUcation by sparsifiable rank-r operators 
has asymptotic complexity 0{rlogr), 

Operation counts: composing two operators. Suppose that M and N aie rank-?* 
operators. Standard muItipHcation of iV and M has complexity O(r^). The complexity 
of injecting N and M into is 0(r^ [log r]^). The action of J^N on J^M has com- 
plexity 0(E„mfc#l-^^v^ • {nY^my^ky) = (n, m, fc)|#| J^Mxr : (ny.mY-ky) = 
(n.m, The second factor is a constant rlogr, while the first when smnmed over 
all nmk is exactly if^\J^N\. Thus the complexity of the nonstandard multiplication 
algorithm, including the conjugation into the basis set a, is 0(#|J^iV| rlogr). Since 
the first factor is r^ in general, the complexity of the exact algorithm is 0(r^ logr) for 
generic matrices, reflecting the extra cost of conjugating into the basis set a. 

For the approximate algorithm, the complexity is 0{#\{J^N)s\ rlogr). For the 
spaxsifiable matrices, this can be reduced by a suitable choice of a to a complexity 
of 0(r^ [log r]^) for the complete algorithm. Since choosing a and cAraluating each 
have this complexity, it is not possible to do any better by this method. 

Reference 

[CW] Appendix IL 

[W2] Appendix V, 
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Appendix IV 
Entropy of a Vector Relative to a Decomposition 

Let V £ H \\v\\ = 1 and assume 
an orthogonal direct sum. We define 

as a measure of distance between t; and the orthogonal decomposition. 

is characterized by the Shannon equation which is a version of Pythagoras* the- 
orem. 
Let 

and Hj give orthogonal decomposition Then 

This is Shannon's equation for entropy (if we interpret as in quantum mechanics 

as the ""probabiUtjr" of t; to be in the subspace J?+). 
This equation enables us to search for a smallest entropy space decomposition of a 
given vector. We need the following H =^ Hi ® H2, 

Hi has two decompositions in JET* or K^, 
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Lemma 1. Let v £ H with \\v\\ and 

Assume also tb&t e^iv{,{ir}) < £^(1^, {K*}) then, if Hi = ® we have 

e^iv, {H'^JJ}) < e\v, {K',L'}) 

Proof, By Shannon's equation 

E^v, {W, L^}) = e^iv, {Hi,H2 }) + ||t;i II V(t;i, {H^}) 
+ ||v2||V(t;^,{27}) 

< e\v, {H^,H2}) + \\vife^{vl{K^}) + |N||V(t;i, {IP}) 
= eHvAK'\L'}). 

Corollary. Assxune £^{vi^{H*}) is the smallest entropy obtsdned for some collection 
of decompositins of Hi and similarly, e^C^g, {1^}) is minimal. Then e^{v; {H^^lJ}) is 
minimal for the direct sum of these collections. 

We consider the following generic example on jD^(R+). 

Let I denote a dyadic interval of the form (2Jfe, 2^k + 1)] , i > 0 fc > 0, and {/« } a 
disjoint cover over (0, oo) consisting of dyadic intervals. We let Hj^ = L^(/a) on which 
we chose an orthonormal basis {e^^^j^.} at fixed (say trig polynomials exp(27ri^)x/^ (as).) 
and consider {e^*^} as an orthonormal basis of L^(R"*"), Thus 

Given a vector v we wish to find 1^ such that {ca,*}) is minimal. In order 
to find lot use a stopping time argument. Starting with intervals of length one 
li = + 1]. We pick a dyadic interval of length two which contedns halves /1J2 of 
length one, Le. J = Ji U J2. We compart 

e^{vxjA4}) with e^ivxjAel'Hei'}) 

and pick the basis given the smallest entropy leading to a cover of L^(R.) by intervals 
of length one and two. We now consider dyadic intervals K of length 4 and compai'e 

€^("XK-{ef}) with e^ivxK-dek}) 


SUBSTITUTE SHEET 


wo 91/18361 


PCr/US91/03S04 


- 44 - 

IV-3 

where form a cover of K by dyadic intervals of length one or two selected previously 
to minimize £^ on each half of if. 

(If the vector function v has bounded support we restrict our attention only to 
dyadic intervals contained in the smallest dyadic interval containing the support of v) 
and continue this procedure up to this largest scale. We cledm that the final partition 
la and corresponding basis provides .the minimal entropy decomposition. In fact, this 
is an immediate consequence of Lemma 1 which shows that given the optimal minimum 
entropy partition any refinement corresponds to worse entropy. 
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V-1 

Higher-Dimensional Best Basis Selection 

Introduction* We introduce a method of coding by orthogonal functions which may 
be used to compress digitized pictures or sequences of pictures, matrices and linear 
operators, and general sampled functions of several variables. The method selects a 
most efficient orthogonal representation fix)m among a large number of possibihties. 
The efficiency functional need only be additive across direct sum decompositions. We 
present a description of the method for pictures, namely functions of two variables, 
using Shannon entropy as the efficiency functional, and mean-square deviation as the 
error criterion. 

Best basis method. In Appendix II b developed a method for generating a Ubraiy 
of orthogonal vectors in R*^ (for large n) together with a notion of admissible subsets 
of these vectors. Admissible subsets form orthonormal bases of wavelet-packets, which 
because of their homogeneous tree structure may be rapidly searched for functiouaJ 
extrema. We can use a family of orthonormal vectors well suited to representhig 
functions of 2 variables. These are products of quadrature mirror filters, as defined 
below: 

Let {pk},{Qk} belong to and define two decimating convolution operators P : 
P ^P,Q:l^ ^l^ 8^ follows: 

OO' oo 

^fk= Pjfj+ik, Qfk= 53 93fi+2k. 

J=— oo 3=1— oo 

P and Q are called quadrature mirror Slters if they satisfy an orthogonality condition: 

PQ* = QP" = 0, 

where P* denotes the adjoint of P, and Q* the adjoint of Q. They are further called 
perfect reconstruction Slters if they satisfy the condition 

P*P + Q*Q=I, 
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where / is the identity operator. These conditions translate to r^trictions on the 
sequences {pk}y{qk}' Let mo, mi be (bounded) functions defined by 

fc=— oo fc=-oo 

Then P, Q are quadrature mirror filters if and only if the matrix below is xmitary for 
aUf: 

This fact is proved in [D], 

Now we can define a number of orthogonal 2-dimensional convolution-decimation 
filters in terms of P and Q. Four of them are simply tensor products of the pair of 
quadrature mirror filters, as in the construction of 2-dimensional wavelets of Meyer 
[M]. 

ij 

P2 ^ Q 0 P, F2V{X, y)^Yl ^(*'i)^2x+iP2y+j 

id 

P3 ^ Q 0 Q» P3v(x, 2^) = 53 J)p2^'^iP^y + i 


There are also pairs of extensions of one dimensional filters: 



Prv{x,y) = 

= 5Zv(iii)*x,iP2y+j = 

i 

QyMl®Q, 

QYv(x,y) 


i 

Px = P<S>I. 

Pxv(x,y) : 


= y)p2x+* 

i 


Qx-oix.y) 


i 


SUBSTITUTE SHEET 


wo 91/18361 PCrAJS91/03S04 

- 47 - 

V-3 

These convolution-decimations have the following adjoints: 

id 

F{v{x,y) = 5Zv(i,j)p2i+xg2j+y 

♦J 


id 3 


The orthogonality relations for this collection are as follows: 

PxPx = PyPy = QxQi = QyQ^y = 
J'-xtQa: = QxPx = ^<^^ = QyPy = O 
/ = P^Px ® = PyPr © QyQy 

By a "picture" we will mean any function 5 = 5{x,y) e /^(Z^). The space 
/^(Z^) of pictures may be decomposed into a partially ordered set W of subspaces 
W(nx,nY,mx>rnY), where mx > 0, my > 0, 0 < < 2"*-^'. and 0 < ny < 2""* . 
These are the images of orthogonal projections composed of products of convolution- 
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decimations. Put W{0, 0, 0, 0) = and define recursively 

W{2nx + i, 2nr + i, mx + 1, my + 1) = F2i^jF2i^jW{nx, ny , m^, my), 
W(2nxynYyTnx + l,my) = PxPxW{nx,nY,mxymY)y 
W{2nx -^l.nY.mx +l,my) - QxQxW{nx,nY,mx,mY), 
W{nx,2nY,rnxyTnY + 1) = PyiVW^(njif,ny,mjic,my), 

W(nx,2ny + l,mx,my + 1) = Q*xQxW{nx,nY,mx,rnY). 

These subspaces may be partially ordered by a relation which we define recursively as 
well. We say W is a precursor of W (write W < W) if they are equal or if W = G*GW 
for a convolution-decimation G in the set {Fq^Fx^F^^Fz^Px^Py^Qx^Qy}- We also 
say that W < W if there is a finite sequence Pi, . . . , of subspaces in W such that 
W ^ ^ • • • X Vii W. This is well defined, since each application of G*G increeises 
at least one of the indices mx or my. 

While {W, is not a tree, it may be made into a tree if we select a subset of 
the relation We will say that W — W(njc,«K» m>:,my) is a principal precursor of 
W' = W{n'j^,n(^,m'x,m'Y) (and write W -< W) if one of the following holds: 

(1) my = mx. and W = G*GW for G € {Fq^Fi.Fz.F^.Px.Py.Qx.Qy}. or 

(2) my < mx. and W = G*GW for G € {Px.Qx}, or 

(3) my > mx, and W = G*GW for G G {Py.Qy}- 

Further, we will say that W -< W if there is a finite sequence Vq, . - . , of subspaces 
in W with W -< Vo -<•••-< -< W\ The relation -< is well defined, since it is a 
subrelation of and it is not hard to see that every subspace W € W has a unique 
first principal precursor. Therefore, {W, -<} forms a (nonhomogeneous) tree, with 
W(0, 0,0,0) at its root. 

Subspaces of a single principal precursor W E W will be called its children. By the 
orthogonality condition, 

(F) W = F^FqW ® F*FiW © F^F2W © F^FsW 

(X) = P^PxW © Q^xQxW 

(Y) = P^PyW © Q^tQyW. 

The right hand side contains all the childien of divided into the groups "F," "X." 
and "Y." Each labelled group of children provides a one-step orthogonal decomposition 
of W. and in general we will have three subsets of the children to choose from. 
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The coordinates with respect to the standard basis of W{nx^nY^ "Zx, my ) form the 
sequence Gi . ^,^^'(0,0, 0,0), where m = max{mjc,my}, and the particular filters 
Gi . . . Gm are determined uniquely by nx and . This is described in Appendix 
n, attached hereto. Therefore we can express in standard coordinates the orthogo- 
nal projections of W(0, 0, 0, 0) onto the complete tree of subspaces W by recursively 
convolving and decimating with the filters. 

Relation with one-dimensional wave packets. Let Wnmk be a one-dimensional 
wave packet at sequency n, scale m and position fc, in the notation of Appensix II. 
Then the element in the {kxyky) position of the subspace W(nx,nv,mx*7nK), at 
the index (x^y), may be written as tfnx,mx,fcx(^)^nr,my.fcv(3/)> which is evidently 
the tensor product of two one-dimensional wave packets. This is easily seen from the 
construction of W(TijK^,ny,mx,Tny): in the x-direction, there will be a total of vix 
convolution-decimations in the order determined by nx^ with the result translated to 
position kxj and similarly in the y-direction. 

We will use the notation w^v for the tensor product of two one-dimensional wave 
packets, with the understanding that the second factor depends on the y-coordinate. 
Since the one-dimensional wave packets are themselves a redundant spanning set, their 
tensor products contain a redundancy of bases for Z^R^. We can search this collection 
of bases efficiently for a best-adapted basis, using any additive measure of information, 
in a manner only slightly more complicated than for the one dimensional case. 

Selecting a best basis. Let S = S{x^ y) be a picture, and let W be a tree of wavelet 
pstckets. Choose an additive measure of information as described in Appendix II, and 
attribute to each node W(nxi'^y,'^x>"^v) the measure of information contained in 
the coordinates of S with respect to the wavelet packets it contains. For example, we 
may use Shannon entropy, 

n{W)= p'logp^ 

where p = {S.Wnxmxhx ® '"^nrmyfe^ >, and W = W(nx,mj^,ny,my ). We will choose 
an arbitrary maximum level in the tree W, and mark all of its nodes as "kept." Pro- 
ceeding up from this level to the root, we will compare 7i{W) for a node W of the tree 

W to the minimum ofJ2w^W'eF^i^')^ 12w^W'ex and Ylw-<W'^Y^i^')' 

If 7i{W) is less, then mark W as ''kepf and mark as ^not kept" all nodes W with 
W -< W; otherwise, markW as "not kept," but attribute to it the minimum of the en- 
tropies of its children. When this procedure terminates at the root, the nodes marked 
"kept" will comprise an orthogonal collection of wavelet packets. 
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It is not necessary to mark all descendants of a '"kept" parent as not kept. The 
complexity of the search algorithm is O(nlogn) if we never change the status of de- 

^ scendants, but instead take for the orthogonal collection only those nodes marked 

"kept" which have no ancestors marked "kept." These may be listed efficiently by 

* indexing the tree in the preorder or depth-first order. 

Error estimates for the best basis. Let H{S) denote the entropy of the picture 
S in the best basis found above. This quantity will be foimd attributed to node 
W(0, 0, 0, 0) at the end of the search. It is related to the classical Shannon entropy Hq 
by the equation 

Wo(5) = ||5||-2W(5) + log||5||2 

The largest exp7io(S) = exp terms of the wavelet packet expansion 

for S contain essentially all the energy of the original picture. Mean square enor 
bounds for specific classes of signals are provided in Appendix IV. 
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CLAIMS: 

1. A method for encoding and decoding an input: signal, 
comprising the steps of: 

applying combinations of dilations and translations of 
a wavelet to the input signal to obtain processed values; 

computing the information costs of the processed 

values ; 

selecting, as encoded signals, an orthogonal group of 
processed values, the selection being dependent on the computed 
information costs; and 

decoding the encoded signals to obtain an output 

signal. 

2. The method as defined by claim 1, wherein said wavelet 
has a plurality of vanishing moments. 

3* The method as defined by claim 1 or 2, further 
comprising transmitting the encoded signals, and receiving the 
transmitted encoded signals before the decoding thereof. 

4. The method as defined by claim 3, further comprising 
storing the encoded signals, and reading the stored encoded 
signals before the decoding thereof. 

5. The method as defined by claim 2, wherein said step of 
applying combinations of dilations and translations of the 
wavelet to the input signal to obtain processed values comprises 
correlating said combinations of dilations and translations of 
the wavelet with the input signal. 

6. The method as defined by claim 2 or 5, wherein 
combinations of dilations and translations of the wavelet are 
designated as wavelet-packets, and wherein the step of applying 
wavelet-packets to the input signal to obtain processed values 
includes: generating a tree of processed values, the tree having 
successive levels obtained by applying to the input signal, for a 
given level, wavelet-packets which are combinations of the 
wavelet-packets applied at a previous level. 

7. The method as defined by claim 6, wherein the steps of 
computing information costs and selecting an orthogonal group of 
processed values includes performing said computing at a number 
of different levels of said tree, and performing said s lecting 
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from among the different levels of the tree. 

8. The method as defined by claim 2 or 7, wherein said 
step of selecting an orthogonal group of processed values 
comprises selecting an orthogonal group having a minimal 
information cost. 

9. The method as defined by claim 8/ wherein said step of 
selecting an orthogonal group of processed values includes 
generating encoded signals which represent said processed values 
in conjunction with their respective locations in said tree. 

10. A method for encoding an input signal, comprising the 
steps of: 

selecting a wavelet having a plurality of vanishing 

moments ; 

applying combinations of dilations and translations of 
the wavelet to the input signal to obtain processed values; and 
selecting, as encoded signals, an orthogonal group of 

processed values. 

11. The method as defined by claim 10, further comprising 

decoding the encoded signals. 

12. The method as defined by claim 10 or 11, further 
comprising transmitting the encoded signals, and receiving the 
transmitted encoded signals before the decoding thereof . 

13. The method as defined by claim 10, wherein said step 
of applying combinations of dilations and translations of the 
wavelet to the input signal to obtain processed values comprises 
correlating said combinations of dilations and translations of 
the wavelet with the input signal. 

14. A method for encoding an input signal, comprising the 

steps of: 

applying combinations of dilations and translations of 
a wavelet to the input signal to obtain processed values; 

computing the information costs of the processed 
values ; and 

selecting, as encoded signals, an orthogonal group of 
processed values, the selection being dependent on the computed 
information costs. 

15. The method as d fined by claim 14, wher in said 
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wavelet has a plurality of vanishing moments. 

16. The method as defined by claim 15, further comprising 
transmitting the encoded signals, and receiving the transmitted 
encoded signals before the decoding thereof. 

17. The method as defined by claim 15, further comprising 
storing the encoded signals* 

18. The method as defined by claim 15, wherein said step 
of applying combinations of dilations and translations of the 
wavelet to the input signal to obtain processed values comprises 
correlating said combinations of dilations and translations of 
the wavelet with the input signal. 

19. The method as defined by claim 15 or 18, wherein 
combinations of dilations of the wavelet are designated as 
wavelet-packets, and wherein the step of applying wavelet-packets 
to the input signal to obtain processed values includes: 
generating a tree of processed values, the tree having successive 
levels obtained by applying to the input signal, for a given 
level, wavelet-packets which are combinations of the 
wavelet-packets applied at a previous level. 

20. The method as defined by claim 19, wherein the steps 
of computing information costs and selecting an orthogonal group 
of processed values includes performing said computing at a 
number of different levels of said tree, and performing said 
selecting from among the different levels of the tree. 

21. The method as defined by claim 15, wherein said step 
of selecting an orthogonal group of processed values comprises 
selecting an orthogonal group having a minimal information cost. 

22. The method as defined by claim 20, wherein said step 
of selecting an orthogonal group of processed values comprises 
selecting an orthogonal group having a minimal information cost. 

23. The method as defined by claim 22, wherein said step 
of selecting an orthogonal group of processed values includes 
generating encoded signals which represent said processed values 
in conjunction with their respective locations in said tree. 

24 . For use in a system which receives an encoded signal 
obtained by: applying combinations of dilations and translations 
of a wavelet having a plurality of vanishing moments to the input 
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signal to obtain processed values; and selecting, as encoded 
signals, an orthogonal group of processed values; a decoding 
method comprising: sequentially applying combinations of 
dilations and translations of a wavelet having a plurality of 
vanishing moments to the encoded signals to decode said encoded 
signals; and outputting the decoded result. 

25. The decoding method as defined by claim 24, wherein the 
encoded signals include identification of the selected orthogonal 
group of processed values, and further comprising the step of 
determining said identification, and performing the sequential 
decoding procedure in accordance with said identification. 

26. Apparatus for encoding an input signal, comprising: 
means for applying combinations of dilations and 

translations of the wavelet to the input signal to obtain 
processed values; 

means for computing the information costs of the 

processed values; 

means for selecting, as encoded signals, an orthogonal 
group of processed values, the selection being dependent on the 
computed information costs. 

27- Apparatus as defined by claim 26, wherein said 
wavelet has a plurality of vanishing moments. 

28. Apparatus as defined by claim 27, wherein said means 
for applying combinations of dilations and translations of the 
wavelet to the input signal to obtain processed values comprises 
means for correlating said combinations of dilations and 
translations of the wavelet with the input signal. 

29. Apparatus for encoding an input signal, comprising: 
means for applying combinations of dilations and 

translations of a wavelet having a plurality of vanishing moments 
to the input signal to obtain processed values; and 

means for selecting, as encoded signals, an orthogonal 
group of processed values. 

30. For use in a system which receives an encoded signal 
obtained by: applying combinations of dilations and translations 
of a wavelet having a plurality of vanishing moments to the input 
signal to obtain processed values; and selecting, as encoded 
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signals, an orthogonal group of processed valu s; ad coding 
apparatus comprising: means for sequentially applying 
combinations of dilations and translations of a wavelet having a 
plurality of vanishing moments to the encoded signals to decode 
said encoded signals; and means for outputtlng the decoded 
result. 
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( FIGURE 17 ) 
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( FIGURE 18 ) 
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